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Abstract 

This work covers methodology of solving QCD evolution equation of the parton 
distribution using Markovian Monte Carlo (MMC) algorithms in a class of models 
ranging from DGLAP to CCFM. One of the purposes of the above MMCs is to 
test the other more sophisticated Monte Carlo programs, the so-called Constrained 
Monte Carlo (CMC) programs, which will be used as a building block in the parton 
shower MC. This is why the mapping of the evolution variables (eikonal variable and 
evolution time) into four-momenta is also defined and tested. The evolution time is 
identified with the rapidity variable of the emitted parton. The presented MMCs are 
tested independently, with ~ 0.1% precision, against the non-MC program APCheb 
especially devised for this purpose. 
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1 Introduction 



The problem of solving numerically the so-called evolution equations of the parton dis- 
tribution functions (PDFs) in quantum Chromodynamics (QCD) is revisited again and 
again in all effort of providing more precise perturbative QCD predictions for the experi- 
ments in the Large Hadron Collider (LHC) and other hadron colliders (e.g. Tevatron). In 
this work we intend to present a methodology of solving QCD evolution equations using 
Monte Carlo techniques for several types of the evolutions, the resulting numerical results, 
including the comparisons with other non-MC numerical methods. 

Two decades ago, when first attempts of solving numerically and precisely the evo- 
lution time dependence of the parton distribution functions (PDFs) according to the 
DGLAP [1] equations were made, it was unthinkable that the Monte Carlo techniques 
could be used for this purpose. It was simply because the computers were too slow by 
several orders of the magnitude. Instead, various faster techniques were developed, based 
mainly on dividing the evolution time into short periods and using discrete grid in x- 
space - they are presently still widely used. Nowdays, with much faster computers, it 
is perfectly feasible to solving numerically the QCD evolution equations with 3-4 digit 
precision for DGLAP and other types of evolutions, albeit it is still much slower than 
with other techniques. 

One may therefore ask the following question: does the MC technique of solving QCD 
evolution equations have some advantages over other techniques which makes it worth to 
pursue in spite of its slowness? In our opinion the MC technique offers certain unique 
advantages. Let us mention the most important ones: Although numerical statistical error 
is usually bigger than for other methods, this error is very stable and robust, not prone 
to any effects related to finite grid or time slicing. Another advantage of the MC method 
is that for many types of partons one may solve the evolution equations for all parton 
types simultaneously, without the need of diagonalizing kernels, that is using PDFs in 
the basis of gluon, singlet quark and several types of the non-singlet quark components, 
and then recombining that back. Finally, the biggest potential advantage is that in the 
MC method one can devise mapping of the evolution time and other variables into four- 
momenta, hence to set-up the starting point for constructing a more realistic treatment 
of the multiparton emission shower, thet is the so-called parton shower MC. Also, the 
extensions from orthodox DGLAP towards more complicated kernels/evolutions featuring 
small x resummations, such as CCFM [2], can be treated with the MC techniques more 
easily than with other methods. 

It should be stressed that this work is closely related with another work of ref . [3] . In 
fact the MC programs of this work are exploited in ref. [3] to test more complicated MC 
techniques of solving evolution equations. The main difference between this work and 
ref. [3] is that here we concentrate on the Markovian class of MC solutions, while ref. [3] 
elaborates on the class of non-Markovian techniques, in which the parton energy fraction 
x and its type / are constrained (predefined). The Markovian MC is better suited for 
the final-state parton cascade while the constrained MC of ref. [3] is better for the initial 
state cascade, for instance in hadron colliders (W/Z boson production). 
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Our paper is organized as follows: In Section 2 we present general form of evolution 
equations and their iterative solutions. In Section 3 we describe in detail three Markovian 
algorithms for solving these equations. Section 4 contains details on evolution kernels and 
form-factors. In Section 5 we give some remarks on Monte Carlo implementations of the 
above algorithms. Section 6 is devoted to the Chebyshev polynomials method of solving 
the evolution equations. In Section 7 we present our numerical results. Finally, Section 8 
summerizes the paper. 



2 General evolution equations 

In this work we shall cover several types of the QCD evolution equations ranging from 
DGLAP [1] to CCFM [2] and their extensions. The generic evolution equation covering 
all types of QCD evolution of our interest reads 



d t D f (t,x) = J2[ duX 



ff/ (t,x,u)D f (t,u). (1) 

The parton distribution function (PDF) is Dj(t, u), with x being the fraction of the hadron 
momentum^ carried by the parton and j being the type (flavour) of the parton. The so 
called evolution time t = \nQ represents in QCD logarithm of the energy scale Q = fi 
determined by hard scattering process probing PDF. The case of the LL DGLAP case [1] 
is recovered with the following identification 

^(«,»,») = ^(^) = ^(^), (*) 

where Pj^,(z) is the lowest order DGLAP kernel. 

In the compact operator (matrix) notation eq. ([1]) reads 

dtD(t) = K(i) D(i). (3) 

Given a known D(t ) at the initial time io, the formal solution at any later time t > i is 
provided by the time ordered exponential 

D(i)=exp(Y K(t')dA D(i ) = G K (i,io)D(io). (4) 
The time- ordered exponential evolution operator reads^l 

/ K(t')dt ; ) =i + X)I[/ dt A>u-MU), (5) 

J t J T.O. n =l i=l J to 



1 Ot, equivalently, the fraction of the eikonal "plus" variable. 

2 Here and in the following we adopt the following conventions J\7=i = A n A n _i . . . A2A1 and 
n"=i J dti = f dt n J dt n -i . . . J dt2 J dt\. The inverse ordering will be similarly denoted with F| i=1 - 
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where (I)/ 2 ji {x 2 , x\) = 8f 2 f 1 8 X2=Xl and the multiplication of the operators is defined as 
follows 

dx' Xf 2f >(t 2 , x 2 ,x')X f/ f 1 (t 1 , x' , Xi) . (6) 

From now on we adopt the following notation^: 

5 x=y = 5(x - y), 9 y<x = l for y < x and 9 y<x = for y > x. 

In the case of the kernel split into two components, K(i) = K j4 (t) + K s (t), the solution 
of eq. (jlj) can be reorganized as follow^ 



B(t) = G KB (t,t )B(t ) + J2 

n= 
oo 

Gn(t, to) = G K s(t, to) + 



71=1 



n / ^ 

n / dt * 



G K s(t, £ n ) 
G K s(t, t n ) 



i=l 

n 



D(t ), 



(?) 

where Gkb is the evolution operator of eq. (EJ) of the evolution with the kernel K B . 
Formal proof of identities in eq. (JTj) can be found in ref. [4]. 

2.1 Resuming virtual corrections 

Monte Carlo method cannot efficiently deal with the non-positive distributions, hence 
resummation of negative virtual part in the evolution kernel is a necessary preparatory 
step. It will be done with help of identity of eq. (ED). We are going resum (negative) 
diagonal virtual part K y = K. B in the kernel 

Xf f >(t,x,u) = X^ f ,(t,x,u) + Xf f/ (t,x,u), Xj f (t,x,u) = -SffS x ^ u Xf f (t,x). (8) 

At this point we do not need to be very specific about Xff,(t, x,u) - we only remark 
that due to infrared (IR) singularity at x = u and / = /' it includes IR cut-off, typically 
u — x > A(x, u, t), causing X v to be also A-dependent. 

Thanks to diagonality of the kernel K y , the corresponding time-ordered exponential 
is easily calculable 

{G K v(M')}//'M = Sff'S^u e-W*\ $ f (t,t'\x) = f dt" X} f (t",x). (9) 
Inserting the above in eq. (JTJ) we obtain 



71=0 



n / 

8=1 Jti-1 



G K v(t, t r 



lK R (U)G K v(ti,t 



i-lj 



8=1 



D(t ). (10) 



More compact notation is obtained with the prescription Yli=k A* = I and Yli=k I dti = \. 
Similarly, we define 9 z<y<x = 9 z<y 9 v<x . 

4 The scope of the index i in Yii ceases at the ciosing bracket, but vaiidity scope of indiced variables, 
like ti, extends until the formula's end. The use of eq. (|6|) is understood accordingly. 
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2.2 Momentum sum rule 



Evolution equations and their time ordered solutions do not require any assumptions about 
the normalization of PDFs and kernels. However, Markovian Monte Carlo methods are 
inherently based on the unitary normalization of the probability distributions (for the 
forward step). Hence, we concentrate on the evolution equations which are supplemented 
with some conservation rule, providing time-independent normalization condition. For 
DGLAP it is the momentum sum rule which is obeyed exactly and is exploited to this 
end (it can also be used for the CCFM class models). It will be also formulated in terms 
of the compact operator formalism. Let us define operator (vector) E acting from the left 
side 



E 



D(t) = / dx^x D f (t,x). (11) 
Jo f 

The momentum sum rule can be stated as the following time conservation law: 

<9 t ED(t) = 0. (12) 
Inserting evolution equation one obtains immediately 

d t E D(t) = EK D(£) = 0. (13) 
The sufficient condition for the above to be true is the following property of the kernel 

EK = 0, (EK) f(u) — ^2 I dx xXf>f(t, x, u) — 0, (14) 

/' Jo 

for any u and /. In particular we have EK V + EK R = 0, from which we can derive 
immediately the virtual part of the kernel 

- (EK v )j(h) = u%} f (t, u) = J2 dxx x ' M ) = (EK R )/(«). (15) 

/' Jo 

From EK = also follows the following usefull identity 

EG K (Mo)=E, (16) 
which provides immediately ED(£) = ED(t )- 



2.3 Markovianization 

The aim is now to transform eq. (jTOj) into a form better suited for the Monte Carlo 
evaluation, using Markovian algorithm. The basic problem is to show how to change the 
integration order from f* Q dt n . . . dt 2 dt x to j* t * dt\ dt 2 . . . J^ 1 dt n , taking into 
account non-commutative character of the product of the kernels in the time ordered 
exponentials. 
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00 r n ft 

D(t) = D(t ) Il t=1 / dt * 9 u>u-i G^iU^K^U) 



It is convenient not only to change the order of the t-integration but also to transpose 
simultaneously (temporarily) both sides of eq. (TTOl) 

G K v(Mn). (17) 

n=0 L Jt ° J 

In the next step we isolate the integration over t±, the outermost one, 
D(*) =D(* )jG K v(Mo) + jT dti G K v(tx,t )K R (tx) 



n=l 



Yl H i=2 l dt i9t i >t i - 1 G K y{U,t i .. 1 )K R (t i 



G K v(t, t n ) 



Closer look into second line in the above equation revealqj that it represents again the 
time ordered evolution operator Gk(£, £1) (with to — > ti). We obtain therefore 



D(t) = D(* ){g k v(Mo) + dh Gkv^^o)^!) G K (t,ti; 
Transposition can be now removed and the integral over t\ is pulled out 



(19) 



D(t)= / dti^G K (t,*i) K fl (t 1 )G K v(t 1) to) + GK^(t,to)^ 1= t J>D(t ) (20) 



to 

The above result can be also presented as an integral equation for the evolution operator 

G K (Mo) = jf dtx {g k (Mi) K*(ti)G K v(Mo) + G K v(*,to)<*ti=t} (21) 

This can be inserted back into eq. (fT9l) many times. The following example shows three 
levels of the nesting 



D(t) = / dtj [ dt 2 

J t V Jt! 



1 2 



dt 3 \ G K (t,t 3 ) K R {t 3 )G K v{t 3 ,t 2 ) + G K v{t,t 2 )S, 



t s =t 



x K R (t 2 )G K v (t 2 , h) + G K v (t, tJS- 



t 2 =t 



xK K (^)G K v (tx,to) + G K v (t, t Q )5 tl =tJ D(t ). 

(22) 

It should be stressed that integration over t\ is now the external one and in the MC it 
will be generated first one. 



After renaming U — > and shifting indices i and n by one. 



If the above nesting is continued to the level N + 1, then one may argue that the 
contribution from the term with Gk(£, £jv+i) for large N decreases like 1/N\, hence in the 
Markovian MC we may use the following formula "truncated" at large fixed N playing a 
role of a dummy technical parameter: 



D(t) = f dtJ f dt 2 [ dtJ 

J t a V Jti L J t 2 I 



*JV-1 



dt N { K R (t N )G- K v(t N ,t N _ 1 ) + G K v(t, tjv-i)^; 



tjv=* 



(23) 



xK fl (t 2 )G K v(t2,*i) + G K v(t,ti)<J t2=t 



xK^^jGKV^.to) + G K y(Mo)5t 1= (jD(to), 

where the integration over t^r+i wa s consumed by S tN+1=t . The above identity will be 
instrumental in constructing MMC algorithm in the following section. 

3 Markovian MC algorithms 

For the Monte Carlo method one needs a (sum of) scalar multi-dimensional integral. For 
the straightforward Markovian algorithm we shall take the following multi-integral 

C = ED(t)=EG K (Mo)D(t ). (24) 

The aim is to generate with the MC method all internal integration variables in the above 
equation. Then, the histogram of the variable x = x n and flavour type / = f n is evaluated 
in the high statistic MC run. Such a histogram is defined by means of inserting Dirac 
delta functions in the above multi-integral: 

jD /( x ) = / dx n dx o (Gk(Mo)) (x n ,Xo) S x=x J ffn D fo (to,x ), (25) 

n=0f n f a J // "' / ° 

where n is the dimensionality of the integral in Gk- 
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(26) 



3.1 Basic formalism 

As a warm-up exercise let us insert D(t) of eq. (1201) into ED(t) and check how the identity 
ED(t) = ED (to) is recovered through explicit integration over t± 

ED(t) = jf dt x {eGk(Mi) K^OGkv^i^oJ+EGkv^^^Jd^o) 

dh | - E K^tOG^fe^o) + EG K v(Mo)^ 1=i jD(t ) 

= | -EG K v(ti,to)|£=ft, + EG K v(t,to)}D(to) = ED (to). 

In the above the most essential was the use of EGk^, ti) = E in the first step, because 
it has allowed to decouple ti-integration from the integrations inside Gk(Mi)- Next, 
EK fl (ti) = — EK. v (ti) was employed, then the evolution equation for Gkv and finally 
G K v(t ,t ) = I was also used. The decoupled inner integrations are explicitly present in 
the following iterative formula 

ED(t) = f dtJ f dt 2 [ dt 3 {... 
J t V Jti L J t 2 I 

. . . jf dt N ^EK R (t N )G K v(t N , tjv-i) + EG K v(t, tjv-i)5tjv=t| 



(27) 



xK R (t 2 )G K v (t 2 , tx) + EG K v (t, tx)5 t2=i 
xK^^OGkv^,^) +EG KV (t,t )5 tl=t )D(t ). 



Again, we would like to stress that the order of the integration starting from t\ and ending 
with tjv is exactly the one which will be realized in the Markovian Monte Carlo algorithm. 

3.2 Straightforward Markovian algorithm 

In the Markovian MC we are going to generate tj, one after another, starting from t\ until 
for certain n, 0<n<N,t = t n+ i is reachecj^l. For this to be feasible in the Markovian 
MC, we have to show with the same algebra as in eq. ( 1271) . that all integrals over tj are 
properly normalized to momentum fraction^ starting with the innermost f* ^ dt^ 

and finishing with outermost J t * dt\. Following the above warm-up example one can show 



6 Maximum number of steps N is large and fixed. Formally, N — > oo is understood. 
7 Unitary normalization is obtained by means of applying 1/xi-i normalization factor. 
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Figure 1: Scheme of the standard Markovian Monte Carlo. 



that the integration over t\ decouples completely from all inner integrations over t2, tjv 
and, therefore, can be generated independently as a first variable in the MC algorithm. 

In the MC generation, whenever 5t n+1 =t term is encountered for the first time, the real 
parton emission chain is terminated. More precisely, for all k > n one may formally define 
tk = t, but they are dummy (not used). 

In ref. [5] it was stated, that every standard (classic) MC algorithm can be reduced 
to a superposition of only three elementary methods: mapping of variables, weighting- 
rejecting and branching. As seen in Fig. [TJ where the above basic MMC algorithm is 
depicted using graphical notation of ref. [5], it is indeed a superposition of branching 



and mapping - every box fi,Xi typically includes more elementary methods (typically 
mappings and branchings). 

What still remains is to define in a more detail the distribution of all three variables 
of ti,fi,Xi of the single Markovian step after generating tj_i, Xi-\ in the preceding 
step: 

t 

1 = — / dU {EK*(iOGKv(Mi-i) + EG K v(Mi-iK=t} 

J { J fi-i 

ti-1 
t 

J dt *{[E, J dxi z i 3C A/*-i(**' a; *' a; i-i) e "* /i - l( *" ti " l) ] +^_i^ =t e-^- 1 ( ^- l) | 



Xi-i 

ti-i 



H— 1 

(28) 

Let us also show the above distribution in a form immediately suitable for the MC gen- 
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eration 



x 



.ft 



(29) 



where virtual form-factor is evaluated using real emission kernels and split into contribu- 
tions from various transition channels according to 



ok). (30) 



t\ El U 

to ^' to o •f' 

Given an uniform random number r G (0, 1), generation of tj is done by means of solving 
the equation r = £/(£«) = e - *^- 1 ^'** -1 ^ for ti, within the range r G [ e ~*/i-i(*'* i - 1 ) ) i]. The 
remaining range r G [0, e - *^- 1 ^'** -1 ^] is mapped into a single point ti = t, that is the point 
where the distribution proportional to 8t=ti resides. Flavour index /j is generated accord- 
ing to normalized discrete probability distribution = d ti ^f 1 f l _ 1 (ti, ti-i) / d ti ^ f i _ 1 (ti, U-i). 
Finally, variable Xi is generated according to the normalized integrand of J dxi in eq. ( 1291 ) . 

The above Markovian MC algorithm of Fig. [1] is completely standard and very well 
known. Practical problem is that the generation of ti, for more complicated kernels 
than in DGLAP case requires numerical evaluation and inversion of the form-factor 
V-i)- Generation of fi is always rather trivial. On the other hand, genera- 
tion of Xi can be also nontrivial. The above problems can be solved, at least partly, by 
more sophisticated versions of the Markovian MC, generally using MC weights, see next 
section. 




Figure 2: Scheme of Markovian Monte Carlo with the internal rejection loop. 



3.3 Weighted Markovian MC algorithms 

In the simplest Markovian MC method with weighted events, which will be referred to as 
an internal loop MMC, the real emission kernel in the distribution used in the generation 
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of Xi is replaced by the simplified one Xj_j-_ %i, — ► K^j. ^(tj, Xj, Xj_i), such that 
3C R < Variables generated according to normalized distribution 

P(Xi) = ^ .} ■ r ^-Xg^fa^Xi-x), (31) 

where 

tl M 

$///(ti, t |u) = y <ft J -^xXf tf (t,x,u) (32) 

<o 

is also simpler than $. The above simplification is corrected by the MC weight 

< = gfelfa^ < 1, (33) 

which is used in the local rejection loop, for every forward step separately, using uniform 
random number r: if r > w\ then generation of Xi is repeated. In this method generation of 
ti is still done using the exact Sudakov form-factor 3>f i _ 1 (t», tj_i|x,;„i). This type of MMC 
algorithm is shown schematically in Fig. [2] and it is essentially a particular realization of 
the basic algorithm of Fig. [TJ In the second method, which will be referred to as a global 
loop MMC the approximate form-factor Qf. ^ti, £j_i|xj_i) = ^hh-xiU-, U-i\ x i-i) is 
used for generation of both ti, /j and Xj. Global correcting weight u> is applied at the 
very end of the Markovian chain. However, the weight is not just n^ 2 ' but it can be 
deduced as follows. According to eq. (|29|) the normalized probability of the forward step 
(i — 1) —> i reads 

dPft ( , , \ r , , 8 



j. -(U-i, fi-i, Xi-i) — — + - 

uXiCtii (34) 

= fl, -^-3C?* ft- x- x- ile - *'*-^**'**- 15 4- & 5 g-^-iCMi-i) 

c<i_i<t,:<t A / i /,_ 1 l , 'i) x i) x «-U t; t vti=tv nfi-i u xi=xi-i e 

Xi—l 

The desired distribution of all variables in MMC event with n emission is 



i=l 



However, in the actual global loop MMC method the distribution of these variables (before 
applying correcting MC weight) is the following 



,-.(n) _ -8 



HtfjL-D-H, (36) 

i=l 

where barring means substitution of exact kernels and form-factors with the approximate 
ones: X —> X, $ — > $. Global correcting MC weight is, therefore, just the usual ratio of 
the exact and approximate distributions 



i i( n ) 

w {n) _ V = g*/ n (*,tn)-* A ,(t,tn) 



J< e */i-i(**.**-0-*/i-i (**•**-!) . (37) 
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Figure 3: Scheme of Markovian Monte Carlo with the global rejection loop. 



The above weight is tested against the random number after the entire MC event 
generation is completed, see the external return loop in Fig. [3j Note that, although 
approximate form-factor and its inverse is used here for generation of ti, the 

exact form-factor is still needed to calculate the global weight. 

Finally, we are going to derive the third method which will be referred to as MMC 
with pseudo-emissions. This method is also known in the literature under the name of 
the Markovian MC algorithm with veto or shortly veto algorithm. In this case we do the 
following modification of the evolution kernel 

Xj fl (t,x,u) = X^ fl (t,x,u) - 5 ff/ 5 x=u X s ff (t,x), 
Xf f ,(t,x,u) = Xf f ,(t,x,u) + 6 ff 6 x=u X s ff (t,x), 

where X^j(t, x) is positive and its magnitude is judiciously chosen as the integral difference 
of the exact kernel X R and the approximate kernel % R > % R (typically the same as in 
the previous methods) 

-i 



(38) 



X b ff (t,u) 



E 



dx - (Xf, f (t, x, u) - Xf, f (t, x, it)) • 



(39) 



In this way we are artificially adding to the real emission kernel finite positive contri- 
butions, which represents real emission of a gluon with exactly zero momentum! This 
extra real emission is compensated immediately and exactly by enlarging negative virtual 
correction. Since the total evolution kernel remains unchanged, 

X fr (t,x,u) = X v ff , (t, x,u) + Xf f , (t,x,u), (40) 

the same time-ordered exponential solution remains valid, D(i) = Gj^(t, to)D(to)- How- 
ever, the difference will occur when resumming virtual negative corrections, because we 
are now resumming the enlarged X v . The basic solution used as a starting point for 
MMC now reads 



m = E 

n=0 

{G^ v {t,t')} ff {x,u 



n 



dti 



_i=l 

5 ff ,5 x=u e-W'W 



® f {t,t'\x) 



D(t c 



dt" X v ff (t",x). 



(41) 



Note that in our older papers describing this method we were denoting Tf = and A/ = <!>/ — <t>y 
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The momentum sum rule still holds and can be used to evaluate modified form-factor 

t\ t\ u 

i |u) = j dt % v ff (t,u) = J dt^ J ~ x Kf, f {t,x,u) 



to t f 

h 



J dt [ j ~ X X /7^' x ' u ) + ^//(^ M ) J ( 42 ) 



to 



J f/ J u 



to 



Obviously, % s was adjusted such that = holds. The immediate important gain is 
that simplified form-factor is used to generate ti, instead of more complicated $/. 



ti=t 



fi —fi-1 









K 



i+l 



t. =t 

i+l 



w <r 



i+P i+l 



fi+lj ' 
%i+l~ *i 



Figure 4: Scheme of Markovian Monte Carlo with pseudo-emissions (veto). 

However, there is one more possible gain from $j = <3>j in the algorithm of generating 
fi and Xi. Due to — > %, the probability of choosing fi should be 

The next should be generated according to X^ ^tj, Xj, including singular part 

proportional to S Xi=Xi l 5ff. However, generating Xi and fi according to this distribution 
can be inconvenient and the following clever trick may be helpful. Let us consider for 
a moment the internal loop MMC algorithm with Pf i = dt i &f i f i _ 1 /dt i <&f i _ l for which Xi 
is generated according to %(xi,...). Give uniform random number r, the fraction of 
MC events obeying r > will be {d ti $f i _ l — d ti & f^) / d ti Now, due to d ti ®f = 
d ti $f this fraction happens to be exactly the same as the fraction of events (d ti ^f i _ 1 — 
®tfi h-i) / ®U® U-\ located in the 8 Xi=Xi _ x 8 fy term! 

One can therefore proceed almost exactly as in the internal loop MMC algorithm, 
that is generate fi according to Pj and Xi according to kernel IK, and next, for events with 
r > wf, instead of repeating generation of fi and Xi for the same ti, one sets /, = and 
Xi = Xi-i (zero momentum real gluon!) and proceeds to generation of the next t i+ i. This 
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completes description and derivation of the algorithm of MMC with pseudo-emissions. 
The advantage of this algorithm is that the numerical evaluation and inversion of the 
possibly complicated exact form-factor $>ff>(t, t'\u) is not required - only the simplified 
version &ff/(t, t'\u) is used. This type of MMC algorithm with pseudo-emissions is shown 
schematically in Fig. HI 

Comparing to other derivations of the veto MMC, in our derivation we reduce veto 
MMC to the standard MMC without the need of repetition of the the explicit resummation 
of the contributions form Gs (which is typically done in the derivations of veto MMC in the 
literature). We believe that the proof presented here is both simpler and more rigorous. 

Finally let us comment on one purely technical point. One may get false impression 
that the above algorithm with pseudo-emissions visualized in Fig. [4] cannot be reduced to 
a superposition of the three elementary methods of ref . [5] . In fact it can be done rather 
easily - the above algorithm is just a variant of the basic algorithm of Fig. [TJ in which the 
branch with W z < r representing emission of another type of real gluon G with exactly 
zero momentum is present. 



4 Kernels and form-factors 

Our main interest is in the CCFM-like evolution with the evolution time being rapidity 
and running coupling constant as dependent on the transverse momentum of the emitted 
gluon. The LL DGLAP will be shown as a reference case, while another with rapidity 
ordering and z-dependent as will be also discussed, as a useful intermediate case between 
CCFM and DGLAP. Running coupling constant 

is taken in the LL approximation. All three types of evolution in this work are essentially 
the same as in ref. [3] , so we shall reduce to a minimum presentation of the corresponding 
three kernels and form-factors. 



4.1 Kinematics 

As already stressed we define explicit mapping of the evolution variables to four-momenta, 
because of possible applications in the parton shower MCs. It will be the same as in ref. [3] 
and is basically that of CCFM model [2]. We define fcf to be the momenta of emitted 
partons, whereas qjf denote the virtual partons along the emission tree. The initial hadron 
carries q£ = 2E^. For each emitted parton we define 

kt = qt 1 -qt = 2E h {xi- 1 -x i )=2E h x i - X {l-z i )\ r ]i = \\n%. (45) 

2 fc ( 

Consequently, the transverse momentum of emitted massless parton reads 

kf = JkfK = Kz~ m = - Zi)2E h e-^. (46) 
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This suggests the convenient definition of the rapidity-based evolution time as 

U = - Vi + ln{2E h ) . (47) 
Now, the transverse momentum of the emitted parton (in units of 1 GeV) becomes: 

kf = e u Xi(l - Zi)/zi = e ti x i - 1 (l - z { ) = e^x^i - x { ). (48) 

4.2 Three types of kernels 

In the following we are going to define matrix elements of the kernels 

(K) ffl (x,u) = X ff/ (t,x,u) = % v ff (t,x,u) + %f f (t,x,u), (49) 

starting with the real emission part Xj^,(t, x, u). It includes implicitly IR cut-off u — x > 
A(x,u). The virtual part %^j,(t,x,u) will be determined unambiguously by imposing 
momentum sum rule. It includes implicitly 5ff5 x=u . We will use as a basic building 
block the real emission part of the LL DGLAP kernel. In order to facilitate numerical 
calculation it is decomposed as follows 

zPf f {z) = S ff [0- z + F ff (z)j + (1 - 6 rf )F f , f (z), (50) 

[z = x/u), with the coefficients Aff and functions Ffif(z) defined in ref. [6]. Let us start 
with pure bremsstrahlung case, real emission part. 

Case (A): DGLAP LL is introduced here as a reference case: 

%f f A \t, x, u) = ^^l l -Pfj{x/u) 9 u . x> _ ue , (51) 

where e is infinitesimally small and z = x/u. 

Case (B): The argument in as is (1 — z)q = (1 — z)e t = k T ju\ as advocated in ref. [7]. 
For the IR cut-off we use A(t,u) = Awe - ': 

%f f B \t^u) = a ^^- x/u)et) l -Pf}{x/u) 6 u _ x>uXe - t . (52) 

Case (C): The coupling constant as depends on the transverse momentum k T = 
(u — x)e t , while for an IR cut-off we choose A(t, u) = A(t) = Ae~*. The kernel reads: 

%ff\t,x,u) = a ^ i{U - x)et) lpP(x/u)e u _ x > Xe - t . (53) 

The generalized kernels beyond the case of the pure bremsstrahlung, for the quark- 
gluon transitions, valid for all three cases X = A, B, C, we define as follows 

xXp?\t,x,u) = 5 ff xXfP(t,x,u) + (1 - y^F^W^AWM, (54) 
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where as in the flavour changing elements have no z- or fc T -dependence and the IR cut-off 
A( x ) is the same as in the bremsstrahlung case. 

Note that the case (C) is fully compatible with the CCFM evolution [2], except that 
for the gluon gluon transitions (bremsstrahlung) the non-Sudakov form-factor assuring 
the compatibility with BFKL [8] is not shown (although it is already present in the MC 
program j§. 

As in ref. [3], for cases (B) and (C), we also introduce slightly modified version of the 
quark-gluon changing kernels elements: 

xXpf\t,x,u) = 5 ff xX*p(t,x,u) + (1 - 5 f/f f siil ~ z)et) F f/f (z)e^ z>Xe - t , 

[ (55) 
xXff \t,x,u) = S rf xXf^(t,x,u) + (1 - 5 rf f s[u{ ^ Z)6) F ff (z)e u _ x>Xe - t , 



with the same arguments of as and IR cut-off as for gluonstrahlung. New variants are 
referred to as cases (B') and (C). One can go back from cases (B') and (C) to (B) and 
(C) by means of applying well behaving MC weight. 

4.3 Form-factors 

Sudakov form-factor resulting from resummation of the virtual part in the kernel was 
defined in eq. The virtual part of the kernel is determined through momentum sum 
rule, see eq. f fl5|) . leading to the following expression 



$/(ti,t |w) = tJ / dt — xXf lf (t,x,u) 
f , Jt Jo u 

= ^2 f dt f —(u-y) Xf, f (t, u - y, u) = V / dt f dz uzXf lf (t, 

r, Jt Jo u v Jt Jo 



(56) 



uz, u) 



where z = x/u and y = u — x = (1 — z)u, see also eq. fl30l) . 
Following decomposition of the LL kernel into three parts 

zPf f {z) = S ff ^ + S ff F ff (z) + (1 - 5 n )F n {z), (57) 

the Sudakov form-factor for practical reasons is split into three corresponding parts: 

$/(Mo|u) = QfthMu) + §)(t x Mv) + ^ c f (h,t \u). (58) 

We show in the following explicit expressions for the above form-factor components 
for most complicated case (C), referring the reader to ref. [3] for simpler cases (A) and 



3 The original CCFM was formulated for pure gluonstrahlung, without quark gluon transitions. 
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(B): 



f h , f 1 , - z)ue l ) A ff „ 

Jto JO 71 L — Z 

2 

= Aff—g 2 (t + \nu,ti + hau]tx), 
Po 



rti ri — z)ue t ) 

J to JO 7r 

^(ti,*o|u) = r dt^^-Y, S dzF f , f (z)9 {1 _ z)u ; 



i>\e~ t i 



where t \ = t i — In A , t\ = t\ — lnA , t\ = In A, while function £ 2 is defined in Appendix 
of ref. [3] in terms of log functions. Two other components $^ and <3>j are evaluated 
numerically for every MC event. This is feasible, provided one integration is performed 
analytically (typically that over v = ln(l — z)) and second integration is done numerically, 
see ref. [3] for the details. 



4.4 Discussion 

In all three cased (A-C) the distributions of the single forward step (parton emission) are 
relatively simple - they are build out of LL DGLAP kernels and as depending on ti Zi 
or Jct- The same distributions enter into form-factor of eq. fl56|) . Practical problems in 
the MC implementations are not so much in the distribution shapes as in the kinematic 
limits. We shall therefore concentrate in the following on this subject. For this purpose 
we will draw the limits of the available phase space in the emission of several gluons 
in the two-dimensional Sudakov logarithmic plane parametrized with variables (k + ,k~) 
and (rj, In k T ) simultaneously. The same integration limits are used in the calculation of 
the form-factors. The translation from evolution times and lightcone variables, U,Xi, to 
rapidities and transverse momenta, (rji, In kf), will be done using mapping of Section. (14.1j) 
in all three cases (A-C)q. 

In in Fig. [5] we start with case (C). The total emission phase space has triangular 
shape and is limited by maximum rapidity (from right) minimum k T (from below) and 
conservation of lightcone plus variable, kf < 2Ef l x^\. Within the above phase space, 
momenta of three emitted gluons k^,i = 1,2,3 are represented by the black numbered 
circles. They are ordered in rapidity. The integration domains for the four consecutive 
form-factors ^^(Ulti-i) in the forward step distributions in eqs. (|34H35|) are also shown 
in Fig. [5] triangle and three trapezoids. 

It is now interesting to compare the phase-space limits in the Sudakov plane between 
the case (C) and the two other cases (A) and (B). The corresponding plots are shown in 
Fig. El The main difference is in the shape of the lower infrared (IR) boundary of the 
emission phase space. In the case (A) of DGLAP it is at the same distance ln(l/e) from 

10 This mapping is primarily adequate for (C). In principle it could be different for (A) and (B). 
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Figure 5: Sudakov plane parametrized two sets of variables (k + , k~) and (?7,lnfc T ). Emission 
of three gluons. Their momenta fcf , i — 1, 2, 3 are marked as black numbered circles. Position 
of the Landau pole marked as dashed red line at k T = A . Phase space limits as in case (C), 
that is CCFM evolution. 



the upper limit, hence rhomboid shapes with the variable widths and constant heights. 
In case (B) the IR limit in k T is lowered by the factor which grows after every 
emission, hence we see the trapezoids with the lower boundary descending deeper and 
deeper into smaller k T . The above illustrates also why the construction of the MMC 
programs evolution type (B) served the role of an intermediate step on the way from 
DGLAP to CCFM. 

Last not least, let us show kinematic limits in the extreme case of one — > 0. This 
limit is treated in CCFM evolution better than in DGLAP, because CCFM in this limit 
coincides with the BFKL evolution [8]. Such a case is illustrated in Fig. [7J where the 
second emitted gluon is very hard, that is with high k T , In fact larger than the scale 
of the hard process. (In this part of the phase space the non-Sudakov form-factor plays 
significant role.) The above kinematic region is properly included in the MMC case (C) 
and also in the CMC of ref. [3]. 
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Figure 6: Sudakov plane parametrized with (k + ,k~) and (77, ln/c T ). Emission of three gluons. 
Their momenta fcf, i = 1,2,3 are marked as black numbered circles. Position of the Landau 
pole marked as dashed red line at (A) Q = A or (B) k T = A . Phase space limits as in case 

(A) and (B) 

5 Monte Carlo implementations 

Studies of the DGLAP evolution, case (A), using the Markovian MCs were already covered 
in refs. [9-11] in particular NLO case was extensively studied in ref. [6]. The main aim 
of these papers was to show that MC method, although slower, is equally precise and 
more versatile as compared to older non-MC techniques, for example grid method based 
QCDnuml6 [12]. These MMCs were also used to test first examples of the constrained 
MCs [13, 14] for DGLAP-type evolution. The main advantage of MC method turns out 
to be very good and stable estimator of the error. The slowness of MMCs is mainly 
the problem in any attempt of fitting deep-inelastic ep data. Here, special pretabulation 
procedures are necessary, see refs. [15]. The above studies of the evolution type (A) using 
MMCs were fairly complete, hence there is no need to repeat them here. 

As already said, we do not show/repeat in this work tests of MMC type (A) and we 
will limit numerical results to comparisons of MMC versus non-MC program APCheb [16] 
for evolutions class (B) and (C). It should be stressed that APCheb was originally working 
only for DGLAP and was upgraded to evolutions type (B) and (C) for the purpose of 
the tests with MMCs. Comparisons of MMC and CMC programs for evolutions type 

(B) and (C) were also done and have been presented in ref. [3]. In this way we have 
in our disposal three completely different programs (sometimes even four) which solve 
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numerically evolution equations of all three types (A), (B) and (C) and provide identical 
results within precision of 0.2%! 

5.1 Reusing MMC type (B) as type (C) 

Historically, the MMC for evolution type (B) with as(e t (l — z)) and IR cutoff 1 — z > 
Ae~* was developed first, before CCFM-like scenario (C). While testing first versions of 
MMC type (C) the following observation was helpful. Examining carefully the propability 
distributions of the single forward step lou^i)^ of eqs. (I34|29l) one may notice that the 

(C) 

whole additional dependence on the variable in u^ ^^ can be absorbed into A and 
Ao: 

A o) = <-W A /z*-i, Ao/x^O. (60) 

Of course, this is the consequence of the relations as(kf) = a^e^l — Zi)Xi-i) and 
kf = e u (l — Zi)Xi_i > A. As a results, we could in the tests of MMC class (C) reuse 
the MMC for a^e^l — z)) by means of reseting A — > X/x^i and A — > A/x,__i, before 
generating each single forward step. The above trick was quite helpful in testing MMC 
class (C), for pure bremsstrahlung. 
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6 Solving evolution equations with Chebyshev poly- 
nomials 

In the previous section the Monte Carlo method for solving the evolution equations was 
presented. For the sake of the comparison, we are going to present an alternative method 
based on the expansion in the Chebyshev polynomials. 

We start from the general form ([1]) of the evolution equations 



dtDf(t, x) — / duK,ff(t,x,u) Df 
f i Jo 



(t,u) (61) 



with the kernel ([8]). The momentum sum rule ( fTTl) imposed on the parton distributions 
allows to determine the virtual part of the kernel (JH1) from the condition ( fl~4"l) . As a result, 
we arrive at the most general form of the evolution equations 

d t (xDf(t,x)) = ^ / duxK,ff,(t,x,u) Djt(t,u) 

r Jo 

- D f (t,x)J2 duu)Cf lf (t,u,x) . (62) 

As an illustration, we consider in detail the evolution equations for the case (C) from 
Section 14.21 The evolution parameter t in this case is related to the rapidity y = 77 of the 
emitted real parton by the relation (T46l . which now reads 

k T = 2E h e- y (u-x) =e t (u-x). (63) 

Here u and x are the longitudinal momentum fraction before and after the emission. In 
the leading logarithmic approximation the real emission kernel takes the form 

lK *, (M , u) = f^P<o>(£) ,(„_,), (64) 

where Pjy; are the leading order splitting functions. In order to avoid the Landau pole 
in a s , we assume that the transverse momenta of the emitted partons are bounded from 
below, 

kr > A > Aqcd ■ (65) 

This further restricts the momentum fractions u for the real emission (see the theta 
function in eq. (1531): 

u > x + Ae~* . (66) 

Changing the integration variable, z = x/u, we obtain for the real emission part of the 
evolution equations ( 1621) 



£ dz a ^ x{1 - Z)/Z) (z) * z D f (t, |) , (67) 
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with the upper integration limit given by 

z R (x,t)= 1 > x. (68) 
1 + Ae 1 

For the virtual part of the evolution equations we interchange u <-» x in the kernel 
( I64j) . Now, the conditions which restrict the u- values read 

u<x, k T — e t (x - u) > A . (69) 

Changing the integration variable, z = u/x, in the second integral of eqs. (1621) . we obtain 
for the virtual term 



D^Zr*™ ^*-*" ^), (70) 



xDf(t, x) y 
where now 

z v (x,t) = 1 - Ae~* > 0. (71) 
In summary, we find the following evolution equations 

z R (x,t) 

rw ^ / „ v-^ f , a s (e t x(l — z)/z) foi , . x ^ / x N 
d t {xD f {t,x)) = J2 dz — p ff( z ) - D f{^- y 

1 X 

z v (x,t) 

-xD f (t,x)Y: f dz as{etx{ l~ z)) zP^jz). (72) 
/' 

These equations are complicated enough to be solved only numerically. In the next section 
we will present the method based on the expansion in the Chebyshev polynomials. 

6.1 Chebyshev polynomial method 

In this method we use the Chebyshev polynomials defined by 

T k (y) = cos(k arccos(y)) , y E [—1,1]. (73) 

The index k = 0,1,2,... denotes the polynomial order. Fixing the order, k = N, we 
consider the equation T N (y) = 0. It has N roots (nodes) given by 

^ = 008^(2-1/2), i = l,2,...,N. (74) 

These roots allow to define the following discrete orthogonality relation for the set of the 
Chebyshev polynomials {T ,Ti, . . . T/v_i}: 

N 

T AVi) T k(Vi) = Cj5 jk , (75) 



i=i 
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where j,k — 0,1, ... ,(N - 1). The coefficients C = N and Cj>i = N/2. 

A function f(x) with x G [a, b] can be approximated with the help of the specified set 
of Chebyshev polynomials in following way 

N 

f(x) « ^VnCnTn-Mx)), (76) 
n=l 

where v\ = 1/2, v n >\ = 1 and y = y(x) is an arbitrary, invertible function which trans- 
forms [a, b] — > [—1,1]. The coefficients c n of the expansion can be calculated from the 
orthogonality relation (1751) . 

2 * 

i=l 

where cCj = y~ 1 (yi) are images of the roots (J74l) in the interval [a,b]. From relations 
(1761) and ( 1771) we see that one only needs the values f(xi) at the Chebyshev nodes to 
reconstruct the function at any other x G [a,b]. This observation is a starting point of 
the method of the solution of the evolution equations (1721) . We simply solve them at the 
Chebyshev nodes x = xy.. 

Therefore, writing eqs. (172|) in a prototype form, 

dD(t,x k ) r*<F*>t) 



dt 



/ dzP(t,z)D{t,x k /z), (78) 



we consider the finite set of the first order differential equations for k — 1, 2, . . . , N. The 
integration on the r.h.s. needs the values of D at any point, thus we use the Chebyshev 
approximation 

2 N N 

D(t,x k /z) n-^^VnDfaxdTn-MT^Mxk/z)). (79) 

n=l i=l 

Substituting into (1751) . we find the following set of equations 

dD{t dt Xk) ^AnQD&Xi) (80) 

i=l 

which can easily be solved numerically [16]. The matrix A k i(f) in these equations, 

2 N rz(x k) t) 
A k i{t) = — y^nTn-iiVi) I dzP(t,z)T n _ 1 (y(x k /z)) , (81) 

n=l 

is computed numerically in the process of finding the solution of eqs. (IHOl) . 

The differential equations which we consider need initial conditions at some initial 
scale D(t = to,x). They are usually specified analytically such that the initial values 
D(t ,x k ) at the Chebyshev nodes are easily calculated. 

The results of the comparison of the solutions of the evolution equations obtained 
using the Monte Carlo and Chebyshev methods are discussed in the next sextion. In 
general, a very good agreement between the results of these two methods is found. 
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Figure 9: For evolution type (C) with a(k T ) plotted are distributions xDf(x) and their ratios 
MMC/APCheb for /=gluon (upper) and / = q + q (lower plot). 
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7 Numerical results 



Although our MMC program was systematically tested against non-MC programs APCheb 
and QCDnuml6 for all evolution types (A-C), we shall show examples of the numerical 
results for the more sophisticated and difficult evolution types (B) and (C). 

Fig. [8] demonstrates distributions xDf(t,x) from MMC and APCheb [16] programs and 
the corresponding ratios MMC/APCheb for the evolution type (B), that is with a((l — z)Q). 
The four curves represent xDf(t, x) for Q — e* = 1, 10, 10 2 , 10 3 GeV. The upper plots are 
for f = G, gluon while lower plots are for / = q + q, quarks and antiquarks taken together. 
The starting quark and gluon distribution at Q = e* = lGeV are defined exactly the same 
as in previous works of refs. [9-11]. Results for all Q = 1, 10, 10 2 , 10 3 GeV were obtained 
in the single MC run of ~ 10 10 MC events. As we see the distributions from two programs 
agree within the statistical MC error of about ~ 0.2%. 

In fig. [9] we show the same type of comparison of MMC and APCheb, but for evolution 
type (C). Again precision agreement within the statistical MC error is reached. 

For the LL DGLAP, case (A), we have reproduced results of ref. [9] with smaller 
statistical errors and removing certain numerical biases which were seen in this paper in 
the gluon case, f = G. We do not show explicitly the corresponding numerical results. 

8 Summary 

We have developed and tested Markovian MC programs for two additional types of the 
QCD evolution equations, in addition to DGLAP. One of them is identical with the 
so-called all-loop CCFM (modulo non-Sudakov form-factor). The corresponding MC pro- 
grams were tested to a high-precision level by means of comparison with the other non-MC 
program APCheb. MMC of this work is also used to test another class of the constrained 
MCs in other independent works, for the same class to QCD evolutions. The aim of these 
exercises is to build basis for the new parton shower implementations. The mapping of 
the evolution variables into four-momenta was also introduced and tested. 
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